import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

import warnings

# Suppress only a specific future warning
warnings.filterwarnings("ignore", message=".*use_inf_as_na option is deprecated.*", category=FutureWarning)
warnings.filterwarnings("ignore", message=".*Pass `(name,)` instead of `name`.*", category=FutureWarning)

data_path = {
    'sentinel': "/kaggle/input/climate-crop/Sentinel_2_NDVI_Timeseries.csv",
    "MOD13A1": '/kaggle/input/climate-crop/NDVI_Timeseries_3km_MyRegion_26july(1).csv',
    "MOD11A1" : '/kaggle/input/climate-crop/MODIS11A1_MyRegion_26july.csv',
    "ERA5_LAND": '/kaggle/input/climate-crop/ERA5_land_3km_MyRegion_26july.csv',
    "CHIRPS":'/kaggle/input/climate-crop/CHIRPS_3km_MyRegion_26july.csv'
}

data_description = {
    "MOD13A1":[
        {
            "feature_name":"NDVI",
            "unit":"",
            "scale":0.0001,
            "offset": 0,
            "min":-2000,
            "max": 10000,
            "description" : "Normalized Difference Vegetation Index"
        },
        
        {
            "feature_name":"EVI",
            "unit":"",
            "scale":0.0001,
            "offset": 0,
            "min":-2000,
            "max": 10000,
            "description":"Enhanced Vegetation Index"
        },
        
        {
            "feature_name":"DetailedQA",
            "unit": "",
            "scale": 1,
            "offset" : 0,
            "description":"VI quality indicators"
        },
        
        {
            "feature_name":"sur_refl_b01",
            "unit": "",
            "scale": 0.0001,
            "offset" : 0,
            "min":0,
            "max": 10000,
            "wavelength":"645mm",
            "description":"Red surface reflectance"
        },
        {
            "feature_name":"sur_refl_b02",
            "unit": "",
            "scale": 0.0001,
            "offset" : 0,
            "min":0,
            "max": 10000,
            "wavelength":"645mm",
            "description": "NIR surface reflectance"
        },
        {
            "feature_name":"sur_refl_b03",
            "unit": "",
            "scale": 0.0001,
            "offset" : 0,
            "min":0,
            "max": 10000,
            "wavelength":"645mm",
            "description":"Blue surface reflectance"
        },
        {
            "feature_name":"sur_refl_b07",
            "unit": "",
            "scale": 0.0001,
            "offset" : 0,
            "min":0,
            "max": 10000,
            "wavelength":"2130nm/2105 - 2155nm",
            "description":"MIR surface reflectance"
        },
        
        {
            "feature_name":"SummaryQA",
            "unit": "",
            "scale": 1,
            "offset" : 0,
            "description": "Quality reliability of VI pixel"
        }
        ],
    
    "MOD11A1":[
        
        {
            "feature_name":"LST_Day_1km",
            "unit": "K",
            "scale": 0.02,
            "offset" : 0,
            "min":7500,
            "max": 65535,
            "description": "Daytime Land Surface Temperature"
        },
        {
            "feature_name":"QC_Day",
            "unit": "",
            "scale": 1,
            "offset" : 0,
            "description": "Daytime LST Quality Indicators"
        },
        {
            "feature_name":"LST_Night_1km",
            "unit": "K",
            "scale": 0.02,
            "offset" : 0,
            "min":7500,
            "max": 65535,
            "description": "Nighttime Land Surface Temperature"
        }   
    ],
    
    "ERA5_LAND":[
        {
            "feature_name":"temperature_2m",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature to which the air, at 2 meters above the surface of the Earth"
        },
        
        {
            "feature_name":"dewpoint_temperature_2m",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature to which the air, at 2 meters above the surface of the Earth, would have to be cooled for saturation to occur. It is a measure of the humidity of the air."
        },
        
        {
            "feature_name":"surface_pressure",
            "unit": "Pa",
            "scale": 1,
            "offset" : 0,
            "description": "Pressure (force per unit area) of the atmosphere on the surface of land, sea and in-land water."
        },
        
        {
            "feature_name":"skin_temperature",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature of the surface of the Earth."
        },
        
        {
            "feature_name":"soil_temperature_level_1",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature of the soil in layer 1 (0 - 7 cm)"
        },
        
        {
            "feature_name":"soil_temperature_level_2",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature of the soil in layer 2 (7-28 cm)"
        },
        
        {
            "feature_name":"soil_temperature_level_3",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature of the soil in layer 2 (28-100 cm)"
        },
        
        {
            "feature_name":"soil_temperature_level_4",
            "unit": "K",
            "scale": 1,
            "offset" : 0,
            "description": "Temperature of the soil in layer 2 (100-289 cm)"
        },

        {
            "feature_name":"volumetric_soil_water_layer_1",
            "unit": "Volume fraction",
            "scale": 1,
            "offset" : 0,
            "description": "Volume of water in soil layer 1 (0 - 7 cm)"
        },
        
        {
            "feature_name":"volumetric_soil_water_layer_2",
            "unit": "Volume fraction",
            "scale": 1,
            "offset" : 0,
            "description": "Volume of water in soil layer 2 (7 - 28 cm)"
        },
        
        {
            "feature_name":"volumetric_soil_water_layer_3",
            "unit": "Volume fraction",
            "scale": 1,
            "offset" : 0,
            "description": "Volume of water in soil layer 3 (28 - 100 cm)"
        },
        
        {
            "feature_name":"volumetric_soil_water_layer_4",
            "unit": "Volume fraction",
            "scale": 1,
            "offset" : 0,
            "description": "Volume of water in soil layer 4 (100 - 289 cm)"
        }
    ],
    "sentinel": [
            {'feature_name': 'QA20',
              'unit': '',
              'scale': 1,
              'offset': 0,
              'description': ''},
            
             {'feature_name': 'B2',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B10',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B11',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B8',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B9',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B7',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'NDVI',
              'unit': '',
              'scale': 1,
              'offset': 0,
              'description': ''},
             
             {'feature_name': 'B3',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             
             {'feature_name': 'B5',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B12',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B4',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B8A',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B1',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
             {'feature_name': 'B6',
              'unit': '',
              'scale': 0.0001,
              'offset': 0,
              'description': ''},
            ],
    
    "CHIRPS":[
        
        {
            "feature_name":"precipitation",
            "unit": "mm/pentad",
            "scale": 1,
            "offset" : 0,
            "min": 0,
            "max": 1072,
            "description": "Precipitation"
        }
        
    ]
}

data_collection = {}

for name,path in data_path.items():
    
    df = pd.read_csv(path)

    for feature_disc in data_description[name]:
        feature_name  = feature_disc['feature_name']
    
        stats = ["_mean","_max","_min","_stdDev"]
        feature_name_stat = [feature_name + s for s in stats ]

        if('scale' in feature_disc):
                scale = feature_disc['scale']
        else:
            scale = 1
            
            
        if('offset' in feature_disc):
            offset = feature_disc['offset']
        else:
            offset = 0

        for fn in feature_name_stat:
            
            df[fn] = df[fn]*scale + offset
        
        if('.geo' in df.columns):
            df.drop(columns = '.geo',inplace = True)

        
        data_collection[name] = df.copy()

####
columns_modis_13A1 = ['date', 'region', 'state','buffer_km','DetailedQA_max', 'DetailedQA_mean', 'DetailedQA_min',
       'DetailedQA_stdDev', 'EVI_max', 'EVI_mean', 'EVI_min', 'EVI_stdDev',
       'NDVI_max', 'NDVI_mean', 'NDVI_min', 'NDVI_stdDev', 'SummaryQA_max',
       'SummaryQA_mean', 'SummaryQA_min', 'SummaryQA_stdDev', 'sur_refl_b01_max', 'sur_refl_b01_mean',
       'sur_refl_b01_min', 'sur_refl_b01_stdDev', 'sur_refl_b02_max',
       'sur_refl_b02_mean', 'sur_refl_b02_min', 'sur_refl_b02_stdDev',
       'sur_refl_b03_max', 'sur_refl_b03_mean', 'sur_refl_b03_min',
       'sur_refl_b03_stdDev', 'sur_refl_b07_max', 'sur_refl_b07_mean',
       'sur_refl_b07_min', 'sur_refl_b07_stdDev']

data_collection['MOD13A1'] = data_collection['MOD13A1'][columns_modis_13A1]
data_collection['MOD13A1']['date'] = pd.to_datetime(data_collection['MOD13A1']['date'])
data_collection['MOD13A1'] = data_collection['MOD13A1'].sort_values(['date','state','region']).reset_index(drop = True)

#####


columns_modis_11A1 = [ 'date', 'region', 'state' , 'buffer_km', 'LST_Day_1km_max', 'LST_Day_1km_mean',
       'LST_Day_1km_min', 'LST_Day_1km_stdDev', 'LST_Night_1km_max',
       'LST_Night_1km_mean', 'LST_Night_1km_min', 'LST_Night_1km_stdDev',
       'QC_Day_max', 'QC_Day_mean', 'QC_Day_min', 'QC_Day_stdDev']

data_collection['MOD11A1'] = data_collection['MOD11A1'][columns_modis_11A1]

data_collection['MOD11A1']['date'] = pd.to_datetime(data_collection['MOD11A1']['date'])
data_collection['MOD11A1'] = data_collection['MOD11A1'].sort_values(['date','state','region']).reset_index(drop = True)


data_collection['CHIRPS'].columns


####
CHIRPS_columns =[ 'date','state','region','buffer_km', 'precipitation_max',
       'precipitation_mean', 'precipitation_min', 'precipitation_stdDev']

data_collection['CHIRPS'] = data_collection['CHIRPS'][CHIRPS_columns]
data_collection['CHIRPS']['date'] = pd.to_datetime(data_collection['CHIRPS']['date'])
data_collection['CHIRPS'] = data_collection['CHIRPS'].sort_values(['date','state','region']).reset_index(drop = True)


####
ERA5_LAND_columns = ['date','state','buffer_km','region','dewpoint_temperature_2m_max',
       'dewpoint_temperature_2m_mean', 'dewpoint_temperature_2m_min',
       'dewpoint_temperature_2m_stdDev', 'skin_temperature_max',
       'skin_temperature_mean', 'skin_temperature_min',
       'skin_temperature_stdDev', 'soil_temperature_level_1_max',
       'soil_temperature_level_1_mean', 'soil_temperature_level_1_min',
       'soil_temperature_level_1_stdDev', 'soil_temperature_level_2_max',
       'soil_temperature_level_2_mean', 'soil_temperature_level_2_min',
       'soil_temperature_level_2_stdDev', 'soil_temperature_level_3_max',
       'soil_temperature_level_3_mean', 'soil_temperature_level_3_min',
       'soil_temperature_level_3_stdDev', 'soil_temperature_level_4_max',
       'soil_temperature_level_4_mean', 'soil_temperature_level_4_min',
       'soil_temperature_level_4_stdDev', 'surface_pressure_max',
       'surface_pressure_mean', 'surface_pressure_min',
       'surface_pressure_stdDev', 'temperature_2m_max', 'temperature_2m_mean',
       'temperature_2m_min', 'temperature_2m_stdDev',
       'volumetric_soil_water_layer_1_max',
       'volumetric_soil_water_layer_1_mean',
       'volumetric_soil_water_layer_1_min',
       'volumetric_soil_water_layer_1_stdDev',
       'volumetric_soil_water_layer_2_max',
       'volumetric_soil_water_layer_2_mean',
       'volumetric_soil_water_layer_2_min',
       'volumetric_soil_water_layer_2_stdDev',
       'volumetric_soil_water_layer_3_max',
       'volumetric_soil_water_layer_3_mean',
       'volumetric_soil_water_layer_3_min',
       'volumetric_soil_water_layer_3_stdDev',
       'volumetric_soil_water_layer_4_max',
       'volumetric_soil_water_layer_4_mean',
       'volumetric_soil_water_layer_4_min',
       'volumetric_soil_water_layer_4_stdDev']



data_collection['ERA5_LAND'] = data_collection['ERA5_LAND'][ERA5_LAND_columns]

data_collection['ERA5_LAND']['date'] = pd.to_datetime(data_collection['ERA5_LAND']['date'])
data_collection['ERA5_LAND'] = data_collection['ERA5_LAND'].sort_values(['date','state','region']).reset_index(drop = True)

# data_collection['sentinel'].columns , ['NDVI_max', 'NDVI_mean', 'NDVI_min', 'NDVI_stdDev']

senitenl_columns =[ 'date','state','region','buffer_km'] +  ['NDVI_max', 'NDVI_mean', 'NDVI_min', 'NDVI_stdDev']

data_collection['sentinel'] = data_collection['sentinel'][senitenl_columns]
data_collection['sentinel']['date'] = pd.to_datetime(data_collection['sentinel']['date'])
data_collection['sentinel'] = data_collection['sentinel'].sort_values(['date','state','region']).reset_index(drop = True)

# !pip install --upgrade seaborn -q
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

df = data_collection['MOD13A1'].copy()
df['month'] = df['date'].map(lambda x: x.month)
df['year'] = df['date'].map(lambda x: x.year)
df = df.sort_values("date")
data_collection['MOD13A1'] = df

df = data_collection['sentinel'].copy()
df['month'] = df['date'].map(lambda x: x.month)
df['year'] = df['date'].map(lambda x: x.year)
df = df.sort_values("date")
data_collection['sentinel'] = df
data_collection['sentinel'] = data_collection['sentinel'].dropna(axis=0)

def ema(data, alpha):
    ema_vals = [data[0]]  # start with the first data point
    for val in data[1:]:
        ema_vals.append(alpha * val + (1 - alpha) * ema_vals[-1])
    return ema_vals

# df_2 = df_1.groupby(['state','region','year','crop' ])[['NDVI_mean','NDVI_max']].max().reset_index()###peak nvdi

##classify crop based on the months
def func_crop_mp(x):
    if(x>=7 and x<11):
        return 'kharif'
    elif(x>=11 or x<=4):
        return 'rabi'
    else:
        return 'nocrop'

def func_crop_MH(x):
    if(x>=6 and x<10):
        return 'kharif'
    elif(x>=10 or x<6):
        return 'rabi'
    else:
        return 'nocrop'

def fun_mark_crop(data):

    state = data.state.unique()
    dfs = []
    
    for s in state:

        if(s =='MP'):
            df =  data[data.state == s].copy()
            df['crop'] = df['month'].map(func_crop_mp)
            dfs.append(df)
            
        elif(s == 'MH'):
            df =  data[data.state == s].copy()
            df['crop'] = df['month'].map(func_crop_MH)
            dfs.append(df)
            
    dfs = pd.concat(dfs,axis =0)
    return dfs
            

data_collection['sentinel'] = fun_mark_crop(data_collection['sentinel'])
data_collection['MOD13A1'] = fun_mark_crop(data_collection['MOD13A1'])  

data_chirps = data_collection['CHIRPS'].copy()

data_chirps['month'] = data_chirps['date'].map(lambda x: x.month)
data_chirps['year'] = data_chirps['date'].map(lambda x: x.year)

# data_chirps['crop'] = data_chirps['month'].map(func_crop)
data_chirps = fun_mark_crop(data_chirps)

data_chirps = data_chirps.sort_values("date")
data_collection['CHIRPS'] = data_chirps

### Merge NDVI and precipitation data

def func(x,df_2):

    offset = pd.DateOffset(months=0, days=15)
    
    x1 = x.iloc[0] 
    x2 = x.iloc[1]

    if(pd.isna(x1)):
        return np.NaN
    elif(pd.isna(x2)):
        x2 = x1 - offset
    else:
        pass

    df_2 = df_2[(df_2['date'] <=x1) & (df_2['date'] > x2) ]

    mean_val = df_2[['precipitation_max','precipitation_mean', 'precipitation_min', 'precipitation_stdDev']].mean()

    mean_val['date'] = x1

    return mean_val
data_ndvi = data_collection['MOD13A1'].copy()

data_ndvi_pre = []

for (s,r),grp in data_ndvi.groupby(['state','region']):
    grp['prev_date'] = grp['date'].shift(1)
    
    df_2 = data_chirps[(data_chirps['state'] == s) & (data_chirps['region'] == r)].copy()
        
    mean_val =  grp[['date','prev_date','state','region']].apply(lambda x: func(x,df_2) ,axis=1)
    
    mean_val['state'] = s
    mean_val['region'] = r

    data_ndvi_pre.append(mean_val)

data_ndvi_pre = pd.concat(data_ndvi_pre,axis = 0)
data_ndvi_precipitation =  data_ndvi.merge(data_ndvi_pre , on =['date','region','state'])

### year on year
columns_to_corr = [ 'NDVI_max', 'NDVI_mean','NDVI_min', 'precipitation_max', 'precipitation_mean','precipitation_min']

data_ndvi_year_crop = data_ndvi_precipitation.groupby(['state','region','year','crop',])[[ 'NDVI_max', 'NDVI_mean','NDVI_min', 'precipitation_max', 'precipitation_mean','precipitation_min']].agg(
    NDVI_peak =pd.NamedAgg(column="NDVI_mean", aggfunc="max"),
    precipitation_mean =pd.NamedAgg(column="precipitation_mean", aggfunc="mean"),
    precipitation_max =pd.NamedAgg(column="precipitation_mean", aggfunc="max"),
    precipitation_std =pd.NamedAgg(column="precipitation_mean", aggfunc="std") )

data_ndvi_year_crop = data_ndvi_year_crop.reset_index()

data_era = data_collection['ERA5_LAND'].copy()
data_era['month'] = data_era['date'].map(lambda x:x.month)
data_era['year'] = data_era['date'].map(lambda x:x.year)


data_era["dewpoint_depression"] = data_era["temperature_2m_mean"]  - data_era["dewpoint_temperature_2m_mean"] 
data_era["surface_air_Temp_Diff"] = data_era["skin_temperature_mean"] - data_era["temperature_2m_mean"]  
data_era["ratio_surface_air_temp"]  = 	(data_era["temperature_2m_mean"]   - data_era["dewpoint_temperature_2m_mean"]) / (data_era["skin_temperature_mean"] - data_era["dewpoint_temperature_2m_mean"])


columns_plot = ['temperature_2m_mean','dewpoint_temperature_2m_mean','skin_temperature_mean','soil_temperature_level_1_mean','volumetric_soil_water_layer_1_mean',
               'dewpoint_depression','surface_air_Temp_Diff','ratio_surface_air_temp'
               ]

data_era_month =  data_era.groupby(['region','month','year','state'])[columns_plot].mean()
data_era_month = data_era_month.reset_index().sort_values(['region','month','year','state'])

data_era_month['year_month'] = data_era_month[ ['year','month']].apply(lambda x: str(x[0]) + ' / ' +str(x[1]) ,axis = 1)
data_era_month = data_era_month.sort_values('year_month')

# data_era_month['diff_temp2m_surface'] = (-data_era_month['temperature_2m_mean'] + data_era_month['skin_temperature_mean'])

# data_era_month

### merge all data sets

def func(x,df_2):

    offset = pd.DateOffset(months=0, days=15)
    
    x1 = x.iloc[0] 
    x2 = x.iloc[1]

    if(pd.isna(x1)):
        return np.NaN
    elif(pd.isna(x2)):
        x2 = x1 - offset
    else:
        pass

    df_2 = df_2[(df_2['date'] <=x1) & (df_2['date'] > x2) ]

    columns = ['temperature_2m_mean','dewpoint_temperature_2m_mean','skin_temperature_mean', 'soil_temperature_level_1_mean', 'soil_temperature_level_2_mean','soil_temperature_level_3_mean',
  'soil_temperature_level_4_mean', 'surface_pressure_mean','volumetric_soil_water_layer_1_mean','volumetric_soil_water_layer_2_mean','volumetric_soil_water_layer_3_mean','volumetric_soil_water_layer_4_mean',
 'dewpoint_depression','surface_air_Temp_Diff','ratio_surface_air_temp']
    
    mean_val = df_2[columns].mean()

    mean_val['date'] = x1

    return mean_val


data_merged = []

for (s,r),grp in data_ndvi_precipitation.groupby(['state','region']):
    grp['prev_date'] = grp['date'].shift(1)
    
    df_2 = data_era[(data_era['state'] == s) & (data_era['region'] == r)].copy()
        
    mean_val =  grp[['date','prev_date','state','region']].apply(lambda x: func(x,df_2) ,axis=1)
    
    mean_val['state'] = s
    mean_val['region'] = r

    data_merged.append(mean_val)

data_merged = pd.concat(data_merged,axis = 0)
# data_nvdi_precipitation =  data_nvdi.merge(data_ndvi_pre , on =['date','region','state'])
data_merged =  data_ndvi_precipitation.merge(data_merged , on =['date','region','state'])

columns = ['NDVI_mean','EVI_mean'] + ['precipitation_mean','temperature_2m_mean','dewpoint_temperature_2m_mean','skin_temperature_mean', 'soil_temperature_level_1_mean', 'soil_temperature_level_2_mean','soil_temperature_level_3_mean',
  'soil_temperature_level_4_mean', 'surface_pressure_mean','volumetric_soil_water_layer_1_mean','volumetric_soil_water_layer_2_mean','volumetric_soil_water_layer_3_mean','volumetric_soil_water_layer_4_mean',
  'dewpoint_depression','surface_air_Temp_Diff','ratio_surface_air_temp']
    
corr = data_merged[data_merged.crop !="nocrop"].groupby(['state','region','crop'])[columns].corr(method = 'spearman')[["EVI_mean","NDVI_mean"]]
corr = corr.loc[[i for i in corr.index if (i[-1] not in ['EVI_mean','NDVI_mean'] ) ]]

columns = ['NDVI_mean','EVI_mean'] + ['precipitation_mean','temperature_2m_mean','dewpoint_temperature_2m_mean','skin_temperature_mean', 'soil_temperature_level_1_mean', 'soil_temperature_level_2_mean','soil_temperature_level_3_mean',
  'soil_temperature_level_4_mean', 'surface_pressure_mean','volumetric_soil_water_layer_1_mean','volumetric_soil_water_layer_2_mean','volumetric_soil_water_layer_3_mean','volumetric_soil_water_layer_4_mean']


columns = ['NDVI_mean','EVI_mean'] + list(selected_feature)


data_merged_ndvi_peak =  data_merged.groupby(['region','state','year','crop']).agg(

    NDVI_peak = pd.NamedAgg(column = 'NDVI_mean',aggfunc = 'max') , 
    precipitation_mean = pd.NamedAgg(column = 'precipitation_mean',aggfunc = 'mean'),
     **{col:pd.NamedAgg(column = col,aggfunc = 'mean') for col in selected_feature}
    
)

data_merged_ndvi_peak = data_merged_ndvi_peak.reset_index()

# data_merged_ndvi_peak[(data_merged_ndvi_peak['year'] < 2024) &
#                                (data_merged_ndvi_peak.state == 'MH')  &
#                                (data_merged_ndvi_peak.region == 'tadklas')&
#                                 (data_merged_ndvi_peak.NDVI_peak > 0.6)].copy()

columns = ['NDVI_peak', 'precipitation_mean'] + list(selected_feature)

columns_best = [col + '_best' for col in columns]
columns_diff = [col + '_diff' for col in columns]

##cross join
data_merged_ndvi_peak_compare =  data_merged_ndvi_peak.merge(data_merged_ndvi_peak,on = ['region','state','crop',],suffixes = ("","_best"),how = 'inner',)
# data_merged_nvdi_peak_compare  = data_merged_nvdi_peak_compare[data_merged_nvdi_peak_compare['NDVI_peak'] < data_merged_nvdi_peak_compare['NDVI_peak_best']]
data_merged_ndvi_peak_compare = data_merged_ndvi_peak_compare.reset_index(drop = True)


data_diff =  -data_merged_ndvi_peak_compare[columns_best].values + data_merged_ndvi_peak_compare[columns].values
df_diff = pd.DataFrame(data_diff,columns = columns_diff )

data_merged_ndvi_peak_compare = pd.concat([data_merged_ndvi_peak_compare  , df_diff] ,axis=1 , ignore_index = False )


###filter out

columns_location = ['region', 'state', 'year', 'crop']

data_merged_ndvi_peak_compare = data_merged_ndvi_peak_compare[(data_merged_ndvi_peak_compare['NDVI_peak_diff'] < -0.1) & 
                                   (data_merged_ndvi_peak_compare['crop'] !='nocrop')][columns_location + ['NDVI_peak','NDVI_peak_best','year_best'] + columns_diff ]

# data_merged_ndvi_peak_compare

# df.groupby( ['region', 'state', 'year', 'crop'])

data_merged_compare = []
data_merged_ndvi_peak_compare = data_merged_ndvi_peak_compare[data_merged_ndvi_peak_compare['NDVI_peak_diff'] <-0.2]

for _,grp in data_merged_ndvi_peak_compare.groupby(['region', 'state', 'year', 'crop']):
    
    grp = grp.sort_values('NDVI_peak_diff',ascending = True)[columns_location +  ['NDVI_peak'] + ['year_best','NDVI_peak_best'] +  columns_diff ]
    grp = grp.iloc[:3] # select top 3 most difference point

    if(grp.shape[0]==0):
        continue
    
    grp = grp.groupby(columns_location + ['NDVI_peak']).agg(
            
            year_comparision = pd.NamedAgg(column = "year_best", aggfunc = lambda x: list(x)),
            NDVI_peak_compare_mean = pd.NamedAgg(column = "NDVI_peak_best", aggfunc = 'mean' ),
        
            **{col : pd.NamedAgg(column = col, aggfunc = 'mean' ) for col in columns_diff }
        
    )
    
    data_merged_compare.append(grp)

data_merged_compare = pd.concat(data_merged_compare,axis=0).reset_index()

pd.set_option("display.max_rows",100)
pd.set_option("display.max_columns",100)

# corr[corr.level_3 == 'precipitation_mean']